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FOREWORD 

This  report  was  prepared  by  the  Metals  Behavior  Branch,  Metals  and 
Ceramics  Division.  This  work  was  initiated  under  Project  No.  2307, 

Task  No.  2307P1 ,  Work  Unit  2307P102,  and  supported,  in  part  under 
Contract  F33615-79-C-5129,  with  Universal  Energy  Systems,  Dayton,  Ohio. 

It  was  administered  under  the  direction  of  Materials  Laboratory,  Air  Force 
Wright  Aeronautical  Laboratories,  Wright-Patterson  Air  Force  Base,  Ohio 
with  Dr.  T.  Nicholas  (MLLN)  as  the  project  engineer.  The  work  was 
performed  by  Dr.  Jalees  Ahmad  of  Battelle  Columbus  Laboratories. 

This  report  describes  work  conducted  from  June  1980  through 
September  1980. 
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SECTION  I 
INTRODUCTION 

In  Solid  Mechanics,  the  most  widely  used  variational  principles  are 
the  virtual  displacement  principle,  the  complementary  energy  principle, 
and  Reissner's  variational  principle  (Reference  1).  In  finite  element 
stress  analysis  the  virtual  displacement  principle  is  most  commonly  used; 
providing  the  usual  displacement  based  elements.  The  complementary 
energy  principle  has  also  been  to  some  extent  exploited  (Reference  2). 

The  third,  and  the  most  general  of  the  three  variational  principles,  is 
Reissner's  variational  principle.  Although  this  variational  principle 
has  remained  largely  unexploited  in  certain  types  of  problems  it  proves 
to  be  quite  useful  (Reference  3). 

The  advantages  gained  by  using  Reissner's  variational  principle  are 
significant  when  the  problems  are  such  that  precise  point  wise  stress 
and  strain  distributions  are  required.  Also  in  problems  where  satisfaction 
of  nodal  force  equilibrium  needs  to  be  repeatedly  checked,  such  as  in  non¬ 
linear  material  problems,  Reissner's  principle  is  advantageous.  The 
most  outstanding  advantage  in  using  Reissner's  principle  is  that  the 
continuity  requirements  on  the  choice  of  shape  functions  are  somewhat 
relaxed  (Reference  4). 

The  disadvantages  associated  with  the  use  of  Reissner's  principle 
need  also  to  be  emphasized.  Firstly,  the  mathematical  foundations  of 
finite  element  formulation  based  on  this  principle  are  not  as  strong 
(References  1-4)  as  in  the  case  of  displacement  based  elements.  Much 
work  regarding  the  proof  and  convergence  of  the  method  still  needs  to 
be  done.  The  second  disadvantage  is  that  the  assembled  coefficient 
matrix  in  this  formulation  is  not  as  well  behaved  as  in  the  case  of 
displacement  formulation  (Reference  5). 

The  purpose  of  this  report  is  to  outline  the  salient  features  of  the 
numerical  formulation  of  finite  element  method  based  on  Reissner's 
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variational  principle.  Numerical  examples,  mainly  of  crack  problems, 
are  presented.  Some  two-  and  three-dimensional  analysis  results  can 
useful  in  fracture  mechanics.  Comparison  with  available  analytical 
results  is  provided. 
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SECTION  II 
FORMULATION 

Reissner's  variational  principle  is  distinguished  by  the  first 
variation  of  the  functional  being  equal  to  zero,  i.e.. 


6  Jr  =  0 


(1) 


where 


JR-  /  -  F,Ui  {Ui,j  *  Uj,i»  d» 

V 

-/  V  Vs  -/  0u  "i  <u,  -  V 


(2) 


in  which 


Wc  is  complementary  energy  function 

a..  is  stress  tensor 

1  J 

U..  are  displacement  components 

U.  .  is  partial  derivative  of  U.  with  respect  to  j 
I  s  J  1 

F.  are  body  forces  per  unit  volume 

T. .*  are  surface  forces  per  unit  area 

v.j  are  direction  cosines  of  the  surface  normal 

U. j*  are  prescribed  displacements  on  Su 

S  is  that  portion  of  the  boundary  on  which  displacements 

u  are  prescribed 

S  is  that  portion  of  the  boundary  on  which  stresses  are 

°  prescribed 

v  is  the  domain. 
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Taking  the  first  variation  of  JR  and  equating  it  to  zero  one  can  obtain 
the  basic  equations  of  elasticity  (Reference  6).  In  simpler  terms  the 
principle  can  be  stated  as  follows: 

JR  =  L'  f  V  (3) 


where 

JR  =  Total  potential  energy 
U  =  Strain  energy 

V  =  Potential  energy  due  to  applied  loads 

In  the  absence  of  initial  strains  and  body  forces 

U  =  f  [o][e]  dv  -  U* 
v 

where,  U*  is  complementary  strain  energy.  Therefore, 

JR  [a][e]  dv  -  U*  +  V  (4) 

Use  of  Equation  1  over  a  subdomain  (element)  leads  to  the  following 
matrix  equation: 


[K][a]  =  [F] 


(5) 


where 


.  -K12^ 

‘  [a] 

'  ^ 

n 

H 

n  x  n 

n  x  m 

[a]  = 

n  x  1 

,  [F]  = 

n  x  1 

[ki2] 

m  x  n 

;  to] 

■  m  x  m 

[r] 

.  m  x  1  - 

[Fr] 
m  x  1 

[Kn]  -f  [Z]T  [D]'1  [Z]  dv  (6) 

v 
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[K12]  =  f  [Z]T  [B]  dv 
v 


f  [l]t  [Y]  ds 

s 

u 


(7) 


[Fr]  =  f  [Y]T  [T]  ds  ,  [Fo]  =  f  [Y]T  [U]  ds 


(8) 


[B ]  is  strain-nodal  displacement  relationship  matrix 

[Z]  is  stress-nodal  stress  relationship  matrix 

[L]  is  boundary  traction-nodal  force  relationship  matrix 

[D]  is  stress-stra i n  constitutive  relationship  matrix 

[Y]  is  boundary  displacement  nodal  displacement  relationship  matrix 

[T]  is  boundary  traction  vector 

[U]  is  boundary  displacement  vector 
[r]  is  nodal  displacement  vector 

[a]  is  nodal  stress  tensor  written  in  vector  form 

The  element  matrix  (Equation  5)  is  then  used  to  form  the  master 
coefficient  matrix  employing  the  usual  finite  element  assembly  procedure. 
A  noticeable  feature  of  Equation  5  is  that  nodal  stress  components  along 
with  the  nodal  displacement  components  appear  as  primary  unknowns  of  the 
probl em. 


Further,  [B]  and  [Z]  matrices  appear  unrelated  and  therefore, 
interpolation  functions  for  displacement  and  stresses  can  be  chosen  quite 
independently.  This  provides  increased  flexibility  to  the  computer  code 
in  terms  of  improving  stress  estimates.  As  for  the  satisfaction  of 
equilibrium  and  compatibility  equations  of  elasticity,  they  are  both 
satisfied  in  an  integral  sense.  This  in  turn  relaxes  the  constant 
derivative  requirement  on  the  shape  functions  (Reference  1). 

By  simple  block  matrix  rules  one  can  write  Equation  5  as  follows: 

-[*„](>]  +  [K12][r]  =  [Fj  (9) 
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and 

[K12]T  [o]  =  [Fr]  (10) 

from  Equation  10  it  is  then  possible  to  recover  the  force  vector  Fr;  a 
procedure  which  is  commonly  necessary  in  problems  involving  nonlinear 
material  behavior. 


For  application  to  problems  in  linear  elastic  fracture  mechanics, 
commonly  sought  parameters  are  the  stress  intensity  factors  (Reference  7). 
The  Mode  I  stress  intensity  factor  (Kj)  can  be  found  by  using  the 
following  equations: 


KI  6  n  .  6  .  36n  KII 

- cos  o  [1  -  sin  -5-  sin  -5—]  -  - — 


sin  j 


[2  +  cos  |  cos  y-]  + 


a  =  —  cos  j  (1  +  sin  |  sin  y)  +  - — —  sin  j 

y  2  2  2  /Zi7  2 


e  39 

cos  2  cos  y  +  .  .  . 


or,  alternatively  by 

U  "  (T  (|i)  cos  f  t1  *  2P  +  sin2  |] 

Kn  /r  \]/2  n  ,  * 

+  q —  sin  |  [2  -  2p  +  cos  j]  +  .  .  .  (13) 

..  KI  /r  \1/2  .  0  r,  ,  2 

v  =  r  (2w)  sin  2  f2  -  2p  - cos  2^ 

KII  /  r  V/2  0  r  t  .  2 

+  G~  (^J  C0S  2  [']  +  2P  +  5in  + 
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where  U  and  V  are  X  and  Y  components  of  displacements  near  the  crack  tip. 
G  is  shear  modulus  and 

p  =  v  =  Poisson's  ratio  for  plane  strain  conditions 

=  i  ■"  for  plane  stress  conditions.  (15) 

The  crack-tip  coordinate  system  is  shown  in  Figure  1.  In  the 
evaluation  of  stress  intensity  factors  Equations  13  and  14  are  commonly 
used.  However,  an  extrapolation  of  values  is  usually  necessary  to  arrive 
at  accurate  results.  Further,  in  three-dimensional  problems  it  becomes 
difficult  to  ascertain  which  value  of  p  should  be  chosen  to  represent  the 
actual  situation. 

In  the  formulation  presented  earlier.  Equations  11  and  12  can  be 
conveniently  used  since  accurate  determination  of  nodal  stress  values 
is  possible. 
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SECTION  III 
NUMERICAL  EXAMPLES 

The  formulation  presented  previously  can  be  used  to  develop  a  variety 
of  finite  element  solution  techniques  by  choosing  various  combinations  of 
displacement  and  stress  interpolation  functions.  Suitable  interpolation 
functions  can  be  chosen  for  particular  applications  and  desired  accuracy 
resul ts . 

For  application  to  problems  in  linear  elastic  fracture  mechanics,  it 
was  determined  through  numerical  experimentation  that  parabolic  shape 
functions  for  displacements  and  linear  interpolation  functions 
(Rference  2)  for  stress  components  provide  accurate  results  both  in 
two-  and  three-dimensional  problems. 

The  numerical  examples  presented  in  the  following  were  solved  by  using 
eight  node  quadrilateral  and  20  node  brick  elements  for  the  two-  and 
three-dimensional  cases  considered.  This  procedure  provided  the  dis¬ 
placement  and  stress  components  at  each  node  from  which  other  quantities 
of  interest  (Kj,  crack  opening  displacement,  etc.)  were  obtained. 

In  the  illustrative  examples  that  follow,  three  different  fracture 
toughness  test  specimen  configurations  are  presented.  An  extensive  study 
of  the  effect  of  element  size  in  the  region  surrounding  the  crack  tip  was 
conducted  and  the  results  indicated  that  by  keeping  the  major  element 
dimension  at  one-fifteenth  of  the  crack  length  the  error  between  the 
computed  Kj  values  and  the  reference  Kj  values  was  consistently  minimum 
(0.1  to  5.0%).  In  all  the  examples  presented  here  the  size  of  the 
elements  surrounding  the  crack  tip  was  therefore  tal.en  to  be  (a/50.0), 
where  "a"  denotes  the  crack  length. 

For  two-dimensional  problems,  Kj  was  evaluated  by  using  the  transverse 
displacement  (V)  of  the  node  closest  to  the  crack  tip  at  6  =  180°  and  by 
employing  Equation  14.  For  three-dimensional  problems,  Equation  12  was 


8 


AFWAL-TR-80-41 82 


used  by  substituting  the  transverse  stress  (a  )  at  0  =  0.0  and 
r  =  (a/50.0).  In  all  the  examples  Poisson's  ratio  of  0.3  was  assumed. 

For  two-dimensional  problems  plane  stress  conditions  were  assumed. 

All  the  results  are  presented  in  nondimensional  form.  Nondimensional 
stress  intensity  factor  (Y),  and  nondimensional  crack  opening 
compliance  (C),  and  nondimensional  load-line  compliance  (6)  are  defined 
separately  for  each  case  considered.  In  the  solutions  British  units 
were  employed  and  specimen  thickness  (B)  was  assumed  to  be  unity. 

1.  TWO-DIMENSIONAL  PROBLEMS 

a.  Three-Point  Bend  Specimen 

For  the  three-point  bend  specimen  shown  in  Figure  3,  a  typical  finite 
element  mesh  is  shown  in  Figure  2.  The  problem  was  solved  for  a/W  ratio 
ranging  from  0.3  to  0.8.  The  results  for  nondimensional  stress  intensity 
factors  (Y)  and  nondimensional  load  point  compliance  (6)  are  shown  in 
Figure  2.  The  referenced  values  were  obtained  from  Reference  8.  Crack 
opening  compliance  values  measured  at  gage  point  (C^)  and  at  the  notch 
mouth  (Cm)  are  included  in  Table  1. 

b.  Compact  Tension  Specimen 

Dimensions  of  the  standard  compact  tension  specimen  are  shown  in 
Figure  4.  The  results  for  a/W  of  0.3  to  0.8  were  obtained  by  the  present 
method. 

The  reference  values  for  Y  and  C(l-l)  were  obtained  from  Reference  8. 
The  load  at  the  pin  hole  was  assumed  to  be  sinusoidally  distributed. 

In  Table  2  the  previous  results  and  compliance  values  C(m-m),  6(p-p), 
and  6(g-g)  at  measured  locations  m,  p,  and  q,  respectively,  are  included. 
This  table  also  shows  the  effect  of  applying  a  concentrated  load  at 
point  p. 
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c.  Double  Notch  Ring  Specimen 

The  geometry  of  a  doubly  notched  ring  under  diametral ly  opposed 
concentrated  tensile  forces  is  shown  in  Figure  5.  This  specimen  can  be 
useful  in  studying  the  crack  behavior  in  thick  shell  type  structures 
usually  encountered  in  nuclear  power  plant,  gun-barrels,  and  in  other 
applications.  Nondimensional  values  for  stress  intensity  factor  and 
compliances  at  measurement  locations  1,  2,  and  at  the  crack  mouth  are 
shown  in  Figure  5.  Table  3  shows  the  actual  finite  element  values. 

2.  THREE-DIMENSIONAL  PROBLEMS 

a.  Three-Point  Bend  Specimen 

The  same  three-point  bend  specimen  of  Figure  3  was  analyzed  using  the 
three-dimensional  formulation  presented  earlier.  Due  to  double  symmetry 
about  X  and  Z  axes  (Figure  6),  only  one-quarter  of  the  specimen  was  used 
in  the  analysis.  Constant  pressure  loading  consistent  with  the  dis¬ 
placement  interpolation  functions  was  applied.  The  problem  was  solved 
for  only  one  a/W  ratio  of  0.375. 

Table  4  gives  the  variation  of  Y,  C  ,  and  <5  across  the  thickness  (B) 
of  the  specimen.  The  definition  of  nondimensional  quantities  is  con¬ 
sistent  with  Figure  3.  The  corresponding  plane  strain  Y  value  for  the 
crack  length  considered  is  7.335  (Reference  8). 

b.  Double  Notch  Ring  Specimen 

A  double  notch  ring  specimen  of  1.0  in  thickness  (B)  and  0.5  in  crack 
length  (A)  was  analyzed  using  the  three-dimensional  formulation.  Inner 
and  outer  radii  of  the  ring  were  chosen  to  be  1.0  in.  and  2.0  in., 
respectively.  Due  to  symmetry  about  X,  Y,  and  Z  axes  only  one-eighth 
of  the  ring  with  appropriate  boundary  conditions  needed  to  be  considered. 
Consistent  constant  pressure  load  was  applied  along  the  inner  radius  at 
right  angles  to  the  crack  line  (Figure  5).  Half  the  thickness  "f  the 
ring  was  divided  into  two  layers  of  brick  elements. 


The  variations  of  Kj,  crack  opening  displacement,  and  load  line 
displacements  across  the  ring  thickness  are  given  in  Table  5. 
Notations  are  consistent  with  Figure  5. 
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SECTION  IV 

SUMMARY  AND  DISCUSSION 

Two-  and  three-dimensional  finite  element  formulations  based  upon 
Reissner's  variational  principle  were  described  and  applied  to  three 
different  fracture  toughness  specimen  geometries.  For  the  three-point 
bend  specimen  and  compact  tension  specimen  the  results  of  Kj  and  load 
line  displacements  agreed  well  with  available  two-dimensional  analytical 
solutions.  In  general,  the  results  for  the  compact  specimen  were  found 
to  be  more  accurate.  One  reason  for  the  three-point  bend  specimen 
results  being  consistently  higher  than  the  corresponding  reference  values 
can  be  the  effect  of  wedge  action  (Reference  9).  This  wedge  action  is 
characterized  by  a  horizontal  force  of  magnitude  (P/n)  being  introduced 
at  the  load  point  when  a  concentrated  force  is  split  in  half  to  solve 
geometrically  symmetrical  problems.  For  further  information  on  wedge 
action,  the  reader  is  referred  to  References  9  and  10.  For  the  double 
notch  ring  specimen  geometry  it  is  observed  that  Kj  varies  linearly  with 
nondimensional  crack  lengths  at  0.3  and  0.8.  This  observation  can  be 
useful  in  designing  test  specimens  for  use  in  parametric  crack  growth 
studi es . 

For  the  three-dimensional  cases  considered,  it  was  found  that  for  a 
three-point  bend  specimen  with  straight  crack  front  the  value  of  Kj  is 
considerably  lower  (23%)  at  the  specimen  surface  than  at  the  middle  plane. 
This  fact  can  be  useful  in  understanding  the  phenomenon  of  cra>  k  front 
"tunneling"  in  thick  specimens.  In  the  case  of  ring  specimen,  using  the 
same  Poisson's  ratio  (=  0.3),  the  variations  in  Kj ,  C,  and  6  follow  the 
same  variation  as  in  the  three-point  bend  specimen. 
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RESULTS  FOR  COMPACT  SPECIMEN  (TWO-DIMENSIONAL) 
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Figure  1.  Crack  Tip  Coordinates 
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A/(  R0-Ri ) 

DOUBLE  NOTCH  RING  TENSION  SPECIMEN 
(Rj/R0)  =  0.5 


Figure  5.  Results  for  Ring  Specimen  (Two-Dimensional) 


